# Gilardi, Fabrizio, "The Temporary Importance of Role Models for Women's Political Representation", American Journal of Political Science
# Code to replicate Figure SI5 (Effect of an additional woman elected in the same municipality on the number of women candidates)
# gilardi@ipz.uzh.ch, 2014-06-23

# Set working directory
setwd("../Data/")

# Load packages
library(Zelig)

# Load data
d <- read.csv("dataset-full.csv")

# Function to compute mode
Mode <- function(x) {
 ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}


# Estimate models

out.3.z <- glm(n.women.cand ~ sl.autob.n.women.elected + n.women.elected.lagged + n.women.cand.previously.elected + n.women.cand.previously.elected + n.men.cand.previously.elected + n.women.cand.lagged + n.people.elected + n.cand + referenda + wom.kr + log.dist.zh + log.pop +  kr.svp + steuerfuss + agglo + land
	+ n.women.elected.lagged*y1978
	+ n.women.elected.lagged*y1982
	+ n.women.elected.lagged*y1986
	+ n.women.elected.lagged*y1990
	+ n.women.elected.lagged*y1994
	+ n.women.elected.lagged*y1998
	+ n.women.elected.lagged*y2002
	+ n.women.elected.lagged*y2006
	+ n.women.elected.lagged*y2010
	, data=d, family="poisson")
summary(out.3.z)

mm.3 <- model.frame(out.3.z)

z.out.3 <- zelig(n.women.cand ~ sl.autob.n.women.elected + n.women.elected.lagged + n.women.cand.previously.elected + n.women.cand.previously.elected + n.men.cand.previously.elected + n.women.cand.lagged + n.people.elected + n.cand + referenda + wom.kr + log.dist.zh + log.pop +  kr.svp + steuerfuss + agglo + land
	+ n.women.elected.lagged*y1978
	+ n.women.elected.lagged*y1982
	+ n.women.elected.lagged*y1986
	+ n.women.elected.lagged*y1990
	+ n.women.elected.lagged*y1994
	+ n.women.elected.lagged*y1998
	+ n.women.elected.lagged*y2002
	+ n.women.elected.lagged*y2006
	+ n.women.elected.lagged*y2010
	, model="poisson", data=mm.3)
summary(z.out.3)


# Vector of election years
yr <- seq(1974, 2010, 4)

# Significance level
sig <- qnorm(0.975)


# Empty matrix to store expected values
pr.l.3 <- pr.h.3 <- pr.d.3 <- matrix(NA, length(yr), 4)
pr.l.3[,1] <- pr.h.3[,1] <- pr.d.3[,1] <- yr


# Empty vetor to store baseline values
baseline.3 <- c()


# Compute expected values fot each election year

for(i in 1:length(yr)){

# Compute baseline value of the spatial lag (modal value)
sl.l <- round(Mode(d$n.women.elected.lagged[d$year==yr[i]]), 2)
baseline.3[i] <- sl.l

# Increase baseline by 0.05 (1 woman out of 20 neighbors)
sl.h <- sl.l + 1

# Vector of explanatory variables with baseline value for spatial lag
x.l <- setx(z.out.3, n.women.elected.lagged =sl.l, sl.autob.n.women.elected =round(Mode(d$sl.autob.n.women.elected[d$year==yr[i]]), 1), n.women.cand.lagged=round(Mode(d$n.women.cand.lagged[d$year==yr[i]]), 1), n.women.cand.previously.elected=round(Mode(d$n.women.cand.previously.elected[d$year==yr[i]]), 1), n.men.cand.previously.elected=round(Mode(d$n.men.cand.previously.elected[d$year==yr[i]]), 1),
	y1978 = unique(d$y1978[d$year==yr[i]]),
	y1982 = unique(d$y1982[d$year==yr[i]]),
	y1986 = unique(d$y1986[d$year==yr[i]]),
	y1990 = unique(d$y1990[d$year==yr[i]]),
	y1994 = unique(d$y1994[d$year==yr[i]]),
	y1998 = unique(d$y1998[d$year==yr[i]]),
	y2002 = unique(d$y2002[d$year==yr[i]]),
	y2006 = unique(d$y2006[d$year==yr[i]]),
	y2010 = unique(d$y2010[d$year==yr[i]])
	)

# Vector of explanatory variables with baseline value + 0.05 for spatial lag
x.h <- setx(z.out.3, n.women.elected.lagged =sl.h, sl.autob.n.women.elected =round(Mode(d$sl.autob.n.women.elected[d$year==yr[i]]), 1), n.women.cand.lagged=round(Mode(d$n.women.cand.lagged[d$year==yr[i]]), 1), n.women.cand.previously.elected=round(Mode(d$n.women.cand.previously.elected[d$year==yr[i]]), 1), n.men.cand.previously.elected=round(Mode(d$n.men.cand.previously.elected[d$year==yr[i]]), 1),
	y1978 = unique(d$y1978[d$year==yr[i]]),
	y1982 = unique(d$y1982[d$year==yr[i]]),
	y1986 = unique(d$y1986[d$year==yr[i]]),
	y1990 = unique(d$y1990[d$year==yr[i]]),
	y1994 = unique(d$y1994[d$year==yr[i]]),
	y1998 = unique(d$y1998[d$year==yr[i]]),
	y2002 = unique(d$y2002[d$year==yr[i]]),
	y2006 = unique(d$y2006[d$year==yr[i]]),
	y2010 = unique(d$y2010[d$year==yr[i]])
	)


# Simulate values for baseline
s.out.l <- sim(z.out.3, x=x.l)
# Simulate values for baseline + 0.05
s.out.h <- sim(z.out.3, x=x.h)
# Simulate difference
s.out.d <- sim(z.out.3, x=x.l, x1=x.h)


# Expected value with baseline
pr.l.3[i,2] <- summary(s.out.l)$stats$"Expected Values: E(Y|X)"[1]
# Lower bound of confidence interval with baseline
pr.l.3[i,3] <- pr.l.3[i,2] - sig*summary(s.out.l)$stats$"Expected Values: E(Y|X)"[2]
# Upper bound of confidence interval with baseline
pr.l.3[i,4] <- pr.l.3[i,2] + sig*summary(s.out.l)$stats$"Expected Values: E(Y|X)"[2]

# Expected value with baseline + 0.05
pr.h.3[i,2] <- summary(s.out.h)$stats$"Expected Values: E(Y|X)"[1]
# Lower bound of confidence interval with baseline 1 0.05
pr.h.3[i,3] <- pr.h.3[i,2] - sig*summary(s.out.h)$stats$"Expected Values: E(Y|X)"[2]
# Upper bound of confidence interval with baseline + 0.05
pr.h.3[i,4] <- pr.h.3[i,2] + sig*summary(s.out.h)$stats$"Expected Values: E(Y|X)"[2]

# First difference
pr.d.3[i,2] <- summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[1]
# Lower bound of confidence interval for first difference
pr.d.3[i,3] <- pr.d.3[i,2] - sig*summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[2]
# Upper bound of confidence interval for first difference
pr.d.3[i,4] <- pr.d.3[i,2] + sig*summary(s.out.d)$stats$"First Differences: E(Y|X1) - E(Y|X)"[2]

}


#pdf(file="fd-poisson-top20-same-f.pdf", paper="special", width=7, height=5.5)
par(mar=c(3.5,3.5,2,1), mgp=c(2.5,0.8,0), cex.axis=1.2, cex.lab=1.2, font.main=1)
plot(yr, pr.d.3[,2], type="n", ylim=range(pr.d.3[,2:4]), xlab="", ylab="Expected increase in nr of women cand.", axes=F)
axis(1, at=yr)
axis(2)
abline(h=0, lty=2)
points(yr, pr.d.3[,2], type="b", cex=1.5)
segments(yr, pr.d.3[,3], yr, pr.d.3[,4])
legend("topright", pch=c(), legend="Effect of an additional woman elected\nin the SAME municipality (t-1)", pt.cex=1.5, cex=1.2, bty="n")
mtext(baseline.3, at=yr, side=1, padj=3, cex=0.9)
mtext("Baseline", at=1974, side=1, padj=3, adj=1.4, cex=0.9)
#points(yr + 0.5, di.2, type="b", cex=1, col="gray")
#segments(yr + 0.5, pr.d.2[,3], yr + 0.5, pr.d.2[,4], col="gray")
#dev.off()



